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ABSTRACT 

I study the formation of Comptonization spectra in spherically symmetric, fast mov- 
ing media in a flat spacetime. I analyze the mathematical character of the moments of 
the transfer equation in the system-frame and describe a numerical method that pro- 
vides fast solutions of the time-independent radiative transfer problem that are accurate 
in both the diffusion and free-streaming regimes. I show that even if the flows are mildly 
relativistic {V ~ 0.1, where V is the electron bulk velocity in units of the speed of light), 
terms that are second-order in V alter the emerging spectrum both quantitatively and 
qualitatively. In particular, terms that are second-order in V produce power-law spec- 
tral tails, which are the dominant feature at high energies, and therefore cannot be 
neglected. I further show that photons from a static source are upscattered by the bulk 
motion of the medium even if the velocity field does not converge. Finally, I discuss 
these results in the context of radial accretion onto and outflows from compact objects. 

Subject headings: plasmas ~ radiation mechanisms: thermal - radiation transfer 



1. INTRODUCTION 

Compton scattering in quasi-spherical accretion flows is thought to be responsible for the X-ray 
spectra of many accreting compact objects. For example, the hard X-ray spectra of weakly-magnetic 
accreting neutron stars can be accounted for if a geometrically thick, hot medium surrounds the 
neutron star and inner accretion disk (see, e.g., Psaltis, Lamb, & Miller 1995). The spectral 
signature of isolated neutron stars accreting from the interstellar medium, which is required for their 
positive identification, is determined mainly by the efficiency of Compton upscattering of thermal 
photons emitted from the stellar surfaces (see, e.g., Zampieri et al. 1995). Moreover, Compton 
scattering plays a major role in the formation of the high-energy spectra of advection-dominated 
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accretion flows, which may be the relevant mode of accretion onto galactic and supermassive black 
holes (see, e.g., Esin et al. 1997). 

Compton scattering in static and moving media, even in plane-parallel or spherical symme- 
try, couples strongly all three independent phase-space coordinates of the photons, namely their 
energy, direction of propagation, and spatial coordinate. For this reason, Monte Carlo methods 
have been used extensively in calculating Comptonization spectra, because of their easy implemen- 
tation and natural treatment of the diverse length scales in an accretion flow (see, e.g., Laurent k. 
Titarchuk 1999 for a recent study; methods and results based on Monte Carlo treatments have been 
reviewed by Pozdnyakov, Sobol, &; Sunyaev 1983 and will not be discussed here). Alternatively, 
Compton scattering can be studied by solving the energy-dependent radiative transfer equation 
or its moments. Implementing such a method is more difficult and often relies on a number of 
approximations for describing the interaction of radiation with matter (see, however, Poutanen &; 
Svensson 1996 and Hsu & Blaes 1998 for complete solutions in static media). However, at the same 
time it is substantially faster and offers significant advantages when the radiative transfer problem 
is coupled to a solution of the fluid dynamics, as well as in treating steep spectra and processes 
that do not conserve photons. 

Over the last twenty five years several authors have derived the photon kinetic or radiative 
transfer equation and their moments for Compton scattering in static and moving media (see, 
e.g., Kompaneets 1957; Babucl-Peyrissac & Rouvillois 1969; Pomraning 1973; Chan & Jones 1975; 
Payne 1980; Blandford & Payne 1981a; Thorne 1981; Fukue, Kato, & Matsumoto 1985; Madej 1989, 
1991; Titarchuk 1994). The moment equations derived by these various sets of authors have been 
widely used in studies of Comptonization by static media (see, e.g., Katz 1976; Shapiro, Lightman, 
& Eardly 1976; Sunyaev & Titarchuk 1980) or by strong shocks and accretion flows onto compact 
objects (see, e.g., Blandford & Payne 1981b; Payne & Blandford 1981; Lyubarskij & Sunyaev 1982; 
Colpi 1988; Riffert 1988; Mastichiadis & Kylafis 1992; Titarchuk & Lyubarskij 1995; Turolla et al. 
1996; Titarchuk, Mastichiadis, & Kylafis 1997), either under the diffusion approximation or with 
closure relations that were specified a priori. Most of the above analyses of Compton scattering in 
moving media were performed in the system frame and to first order in the electron bulk velocity. 
On the other hand, Schmid-Burgk (1978), Zane et al. (1996), and Titarchuk & Zannias (1998) have 
solved the general-relativistic radiative transfer equation for Compton scattering without the need 
for specifying a priori closure relations. 

In Paper I (Psaltis & Lamb 1997) of this series, the time-dependent photon kinetic and radiative 
transfer equations were derived, as well as their zeroth and first moments that describe absorption, 
emission, and spontaneous and induced Compton scattering in static and moving media, correcting 
various errors in the literature. The system-frame equations that were derived are valid to first 
order in e/mg and Te/me, and to second order in V , where rUe and Tg are the electron mass and 
temperature, e is the photon energy, and V is the bulk velocity of the electrons in units of the 
speed of light; the fluid-frame equations that were derived are valid for arbitrary values of the bulk 
velocity V . Using these equations it was argued in Paper I that the effects of Comptonization by the 
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bulk electron velocity that are described by the terms that are second-order in V are qualitatively 
different than and can become at least as important as the effects described by the terms that are 
first-order in V, even when V is small. As a result, these second-order terms should generally be 
retained (see also Yin &: Miller 1995). 

In this paper, I describe a numerical algorithm for the solution of time-indcpcndcnt radiative 
transfer problems in systems with spherical symmetry. Even though the natural reference frame 
for solving transport problems is the one comoving with the flow, I work in the system frame, 
which can be any inertial frame in which the central object is at rest. This is important for the fast 
algorithm described below, because the moment equations in the system frame are always parabolic 
(as shown in §3) in contrast to the fluid-frame equations (Turolla, Zampieri, & Nobili 1995; Smit, 
Cernohorsky, &: Dullemond 1997; see also Korner & Janka 1992). Moreover, in problems of mass 
accretion onto a compact star, the boundary conditions are often more easily implemented in the 
system frame. I also neglect any general relativistic effects, which can be shown to introduce 
only small quantitative corrections to the emerging spectra (Papathanassiou &; Psaltis 2000), and 
truncate the transfer equation keeping only terms up to second order in V. 

The algorithm used is a generalization of the method of variable Eddington factors (Mihalas 
1980; see also Mihalas 1978, p. 157-158). It is therefore accurate for systems of arbitrary optical 
depth and does not require a priori specification of closure relations. It is based on the iterative 
solution of both the transfer equation and the systems of its first two moments and hence the 
validity and accuracy of the solution can be verified explicitly at the end of each run. The method 
also requires only a small number of iterations in order to converge to the solution and is therefore 
ideal for future extensions to time-dependent and multi-dimensional transport problems. 

In §2, I summarize the results from Paper I and introduce the notation. In §3, I discuss the 
mathematical character of the system of equations and in §4, I present the numerical algorithm. In 
§5, I describe the results of numerical solutions of the transfer equations for a variety of situation. 
Finally, in §6, I briefly summarize the implications of the above for Comptonization spectra in 
static and moving media. 

2. THE ELECTRON GAS AND RADIATION FIELD 

Units — Throughout this paper I set h = k-Q = c = 1, where h is Planck's constant, fee is 
Boltzmann's constant, and c is the speed of light. I also normalize all the spatial coordinates to the 
inner radius of the flow, the electron temperature and photon energy e to the electron rest mass 
me, and the electron density to its value at the inner radius of the flow. Finally, I normalize 
the emission and absorption coefficients r]{r, e) and x(r, e) to the inverse of the electron scattering 
mean-free path in the Thomson limit. 

The electron gas — Let u be the system-frame three-velocity of a given electron. I define the 
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local bulk velocity V of the electrons measured in the system frame as 

V={u), 



(1) 



where the sharp brackets indicate the average over the local electron velocity distribution. I also 
assume that the electron velocity distribution in the fluid frame is a relativistic Maxwellian and 
therefore 

(n^) ~ y2 ^ 3Te , (2) 

where I have neglected terms of order V'^T^ and higher. 

The radiation field — I describe the radiation field at any position in the system frame with 
coordinate vector r using the monochromatic specific intensity I{r,l,€), where I is the direction of 
propagation and e is the photon energy. Here I consider only unpolarized radiation and hence I 
have suppressed the dependence of the specific intensity on polarization. Because of the assumed 
spherical symmetry of the problem, the radiation field depends only on the radial distance from 
the center of symmetry (i.e., the r-component of the vector r), on the angle 6 between the radial 
direction and the direction of propagation (i.e., cos^ = / ■ f/r), and on the photon energy e. 

I define the first five moments of the monochromatic specific intensity as 
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The non-zero components of the first five moments of the specific intensity for systems with spherical 
symmetry are given in Appendix A. 

The transfer equation and its moments — Following Mihalas (1980; see also Mihalas 1978), I 
solve the radiative transfer equation in a plane that contains the center of symmetry of the system, 
between the inner and outer radii rin and rout of the flow, using the coordinates z = r cos 9 and 
p = r sinO. 



The radiative transfer equation along a ray of constant impact parameter p is 

± I^(z,p,e) = (1 -2e)I±(z,p,e) - 5±(z,p,e) , 



(8) 



where the signs '+' and '— ' correspond to the equations for the outgoing and incoming rays re- 
spectively and I have neglected the effects of induced Compton scattering. In equation (8) I have 
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used the optical depth along the ray defined as dTz = —UeCrTi^ + xi''^^ ZT^)]dz, where = at 
the outer boundary and (Tt is the angle-integrated Thomson scattering cross section. {z,p,e) is 
the generalized source function for absorption, emission, and Compton scattering in systems with 
spherical symmetry and is given in Appendix B. 

The equations for the moments of the specific intensity can be derived either directly from 
equation (8) or from equations (34) and (40) of Paper I. Written in compact form and using the 
dimensionless quantities defined above, the moment equations are^ 

AiJ + A2€d,J + A3e^d^^J + AiW + A5€d,H'^ + dr^H'- = Ci (9) 
BiJ + B2ed,J + rdrj + B^H'' + Bied,H'' + B^e^dlH' = C2 , (10) 

where = d/de^ = (9^/3e^, dr^ = d/drr^ and = — cJT?i-e(nn)(?' ~ ^'out)- When the electron 
density is uniform, then the dimensionless quantity Tr is equal to the electron scattering optical 
depth in the radial direction measured from rout- The coefficients in equations (9) and (10) are 
given in Appendix C. 

Boundary Conditions — The boundary conditions for both the transfer equation and the system 
of its moments depend on the specific problem under study. For the model problems discussed in 
§5, I shall assume that the flow is not illuminated from the outside, i.e., 

/-(rout,^,e) = (11) 

at the outer boundary, and specify the radiation flux at the inner boundary as 

/+(rin, 9, e) = 7-(rin, 9, e) + AHl{e) . (12) 

For solving the moments of the radiative transfer equation I shall specify, for all photon energies 
e, the flrst moment of the monochromatic specific intensity at the inner radius, i.e., 

i?'-(n„,6) = i75(6) (13) 

as well as the ratio of the first to the zeroth moment of the monochromatic specific intensity at the 
outer radius 

The quantity fc(e) is not known a priori but is calculated in a self-consistent way in the iterative 
method described below. I shall also set at all radii 



e-*0 

and 



lim J(r, e) = lim J(r, e) = (15) 



lim iJVr, e) = lim H'^ir, e) = , (16) 
which reflect the fact that the photon occupation number and the radiation flux are flnite quantities. 



^Note that the B3 term reported in Paper I (Psaltis & Lamb 1997) contained a typo that has been corrected here. 
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3. MATHEMATICAL CHARACTER OF THE MOMENT EQUATIONS 

In the method of variable Eddington factors, the system of partial differential equations (9) 
and (10) is solved iteratively with the radiative transfer equation (8), given a set of Eddington 
factors computed from the previous iteration. The boundary conditions and the method used 
for the solution of the system depends on the character of the partial differential operator. For 
example, the system of moment equations in the frame comoving with the flow (or the local Lorcntz 
frame when gravitational effects are taken into account) may be hyperbolic or elliptic, depending 
on the velocity field (Turolla et al. 1995). In this section I follow the procedure outlined by Ames 
(1992, pp. 8-12) to show that the system-frame, energy dependent moment equations always form 
a parabolic system of partial differential equations, thus simplifying the implementation of the 
numerical algorithm. 

I define the quantities L = d^J and = d^H^, which together with the partial differential 
equations (9) and (10), and the equations for the differentials dJ, dH^, dL, and dAf, form an 8 x 8 
system of algebraic equations for the eight derivatives 
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Ci - AiJ - A2eL - A^H'' - A^eW 
(72 - BiJ - B2eL - BsH - B^eM' 

dJ-Lde ■ ^ ' 

dH"- - Wde 
dL 
dM^ 

Equating the determinant of the coefficient matrix to zero, I obtain the characteristic equation for 
the system of equations, 

AgSse^ (dr,)^ = . (18) 

The coefficients A^ and B5 are proportional to the electron temperature and the square of the bulk 
velocity (see Appendix C) and hence for the problems under consideration they are never zero. 

Therefore, for a finite electron energy e, the system of partial differential equations (9) and (10) has 
only one multiple, real characteristic direction, i.e., dTr = 0, and hence is always parabolic. The 



- 7- 



difference with respect to the moment equations in the comoving frame arises from the presence of 
second order energy derivatives in the differential operator. 



4. NUMERICAL METHOD 

Equation (8) is an integro-differential equation for the specific intensity of the radiation field 
and therefore any numerical solution of this equation is not trivial to obtain. On the other hand, 
the system of equations (9) and (10) depends on two variable Eddington factors and an unknown 
outer boundary condition, which can be computed only when the solution of equation (8) for the 
specific intensity is obtained. Here, I employ a numerical algorithm that is a generalization of the 
method of variable Eddington factors (Mihalas 1980; see also Mihalas 1978, pp. 157-158) in order 
to solve iteratively equation (8) for the specific intensity of the radiation field and the system of 
partial differential equations (9) and (10) for its zeroth and first moments. This algorithm has been 
proven to be very efficient and gives a solution to the complete transfer equation that is limited 
only by numerical accuracy, independent of the optical depth of the medium. 

In this method, I first solve the moments of the transfer equation using an initial guess for the 
variable Eddington factors and for the outer boundary condition. I then compute the generalized 
source function using the calculated moments of the specific intensity, solve the radiative transfer 
equation, and use this solution to update the Eddington factors. I repeat the whole procedure until 
convergence is achieved. 



4.1. Solution of the Moments of the Transfer Equation 

The moments of the transfer equation depend only on radius and photon energy, when the 
variable Eddington factors are known. Therefore, I discretize the domain of solution and the 
differential operators of equations (9) and (10) on a two-dimensional mesh of Nt grid points over 
the variable and Neu grid points over the photon energy e. 

The mesh of discrete grid points over the variable must resolve the rapid change of the 
variable Eddington factors with optical depth near both boundaries of the domain of solution. In 
order to achieve this, I define the mesh points to be equidistant in the quantity 

{WT+{l-W)logT, T<\Tm , . 

\W{Tm-T) + {l-W)\og{T^-T) T > ' ^ ' 

where r is the total optical depth (and not the quantity r,.) in the radial direction measured from 
^out) '^m = T{fin) is the maximum radial optical depth in the medium, and is a parameter 
that allows a combination of a logarithmic grid near the boundaries with a linear grid in the 
interior. Because the spectra that emerge from Comptonizing media are usually power laws, I use 
a logarithmic mesh of discrete grid points over photon energy e. 
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The system of moment equations involve differentiation over two density-like quantities, the 
radiation energy density J and the first Eddington factor f"', and two flux- like quantities, the 
radiation energy flux H'^' and the second Eddington factor g"' . I discrctize density- like quantities 
in a shell-centered fashion and flux-like quantities on the grid points. Let Fij be the value of a 
physical quantity at the i-th grid point in the variable and the j-th grid point in photon energy. 
Away from the boundaries, I differentiate density-like quantities over the variable on the grid 
points as 

—\ = 2Ei+lI2A.lIhllM. (20) 

dTrJi^j {Tr)i+1 - {Tr)i-1 



and flux-like quantities as 
in an explicit or fully implicit scheme. 



In both equations I approximate the first derivative of a physical quantity F with respect to 
photon energy by 



and the second derivative with respect to photon energy by 



d'^F \ ^ 2 f Fjj+i - Fjj _ Fjj - Fjj^i > 



where the index 'z' in the quantity Fij may also be fractional, corresponding to a shell center, as 
needed by an implicit differencing scheme over the variable r^. The energy differencing scheme (22)- 
(23) cannot be used directly in differencing the radiation energy density and flux in equations (9) 
and (10) because it does not guarantee particle conservation when only scattering is taken into 
account (Chang & Cooper 1970). This is important when the energy spectrum has reached quasi- 
equlibrium and small numerical errors lead to significant deviation from particle conservation (see 
eq. [10] in Chang & Cooper 1970). A better differencing scheme can be deduced by examining the 
properties of the differential operator at that limit. 

The solution of the moment equations reaches quasi-equlibrium when the photon mean free 
path is very small compared to any characteristic length scale in the system. In that limit, = 
1/3 -|- 0{V), and when only processes that conserve photon number are operating, the moment 
equation of zeroth order becomes 



Or -\ = UeCTTOe 

r I e 



VH' + (e - 3Te - V')J + ( Te + — ) ed,J 



(24) 



where Ug is the dimensional electron density. The photon number flux is simply equal to = H^/e 
and integrating equation (24) over photon energy gives 
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i.e., photon conservation. Following Chang k, Cooper (1970), I define a generalized photon current 
J^{r,e) as 

J^{r, e) = VH"- + (e - ST^ - V^)J + (j^e + ed,J (26) 

and write the zeroth and first moment equations as 

A'^J + A'2ed,J + A'^e^d'^,J + A'^H'' + dr,H'' = Ci+need^J^ (27) 
BiJ + B2edeJ + rdrrJ + B3H' + B^edeH' + B5e^d^H' = C2 , (28) 

where the coefficient are given in Appendix D. Differencing equations (27)-(28) according to the 
scheme (20)-(23) results in strict photon number conservation in the limit of small photon mean 
free paths. 

A final point with the energy differencing is related to the requirement that the photon energy 
density is everywhere positive. This introduces a restriction on the spacing of the energy grid, 
depending on the differencing scheme. For the scheme used above it can be shown that the re- 
quirement for a positive photon energy density results in the restriction (cf. eq. [24] and Chang &; 
Cooper 1970, eq. [14]) 

Ae /1\ Te + yV3 ^fl/6, e«3re + y2 

e ^ V2y |e-3(re + yV3)p I (^)(re + yV3) , €»3re + y2 ' y ) 

which can be easily met. 

I solve the system of equations (27) and (28) for the zeroth and first moments of the specific 
intensity of the radiation field, J and , in Nj- grid points in optical depth and A^En grid points 
in photon energy. Therefore the system of equations has Nr x Aeh unknowns. In all the interior 
grid points (1 < i < A^; 1 < j < Ae^) I solve both equations (27) and (28) using the differencing 
scheme of equations (20)-(23). For the inner spatial boundary, i.e., for i = 1 and for all j, I solve 
equations (13) and (28), whereas for the outer spatial boundary, i.e., for i = and for all j, I solve 
equations (14) and (27). Finally, at both energy boundaries, i.e., for all i and for j = l,NEn, I set 
both the zeroth and first moments of the specific intensity equal to zero, according to boundary 
conditions (15) and (16). 

The resulting system of equations is linear in both the zeroth and first moments of the specific 
intensity. Defining the quantity Pi as 

Pi = Ji,i ; P2 = Hli ; ... P2N^-i = JNr,i P2Nr = Hn^,i etc. (30) 

the system of equations takes the form shown in Figure 1. The matrix of coefficients is a sparse 
matrix and can be solved efficiently with the biconjugate gradient method (see, e.g.. Press et al. 
1992, pg. 209) with a preconditioner matrix. I choose the preconditioner matrix to be the coefficient 
matrix, with the elements at the outer two diagonal bands set equal to zero. The preconditioner 
matrix can then be inverted efficiently with the LU-decomposition method (see, e.g.. Press et al. 
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Fig. 1. — Schematic representation of the Unear system of the discretized zeroth and first moments of the radiative 
transfer equation. 



1992, pg. 202). In this method, in order to avoid round-off errors, I scale aU the rows in the 
coefficient matrix so that the highest value of the elements in each row is unity. Depending on the 
complexity of the particular problem, the system of difference equations can be solved to a fractional 
accuracy of 10~^ within less than fifty iterations, which can be performed in a few seconds on a 
current workstation. 



4.2. Solution of the Radiative Transfer Equation 

In the iterative method of the variable Eddington factors, the transfer equations (8) is solved 
at every grid point in the {p — z) plane for the boundary conditions (11) and (12). Because the 
generalized source function is calculated using the results from the previous iterations, the solution 
of equation (8) is just the formal solution that can be obtained analytically, i.e., 

I-{T,,p,e) = r S{t,p,e)e^'-^'^('--^Ut 
Jo 

I+iT,,p,e) = I+(^min,P,e)e(^-2^)(^^-^^-^") 

fTzmin 

+ 5(t,p,e)e(i-2^)("^-*)di, (31) 
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Fig. 2. — Maximum fractional error in the radiation energy density after each iteration, for the test problems 
discussed in §4.3. 



where 



7"+fr • T) f ^ - J ^ (^min,P, e), P > Hn /or,^. 



and Zmin = or Zmin = {i^in ~ P^)^^"^ when p > nn or p < Tin respectively. In calculating the 
integrals in equation (31), I use a piece- wise linear interpolation of the source function through 
each shell and evaluate the energy derivatives according to equations (22)-(23). 



4.3. Convergence and Validation 

The rapid convergence of the method of variable Eddington factors is preserved in the current 
implementation for the problem of Compton scattering. This is shown in Figure 2, where the 
maximum fractional error in the radiation energy density is plotted for each iteration, for two 
uniform media with different scattering optical depths. For the model problems shown, the outer 
radius of the scattering medium is set to 3r[^, its electron temperature to 10 keV, and the energy- 
dependence of the radiation flux at the inner boundary to be that of a blackbody with a temperature 
of 0.5 keV. The converge is typically faster for problems with higher optical depths, because I have 
used f^^ = 1/3 as an initial value for the first variable Eddington factor, which is closer to the real 
value at the limit of high optical depth. 
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The numerical method described in this section also allows for a consistency check between the 
solution of the transfer equation and the solution of the moment equations that can be performed 
for each individual problem. For this test, equation (8) can be solved iteratively without using the 
solution of the system of equations (9)-(10), by starting with an assumed generalized source function 
[e.g., S{z,p, e) =0], solving the transfer equation for the specific intensity, updating the generalized 
source function, and repeating until convergence is achieved. The moments of the specific intensity 
can then be obtained in this way and compared directly with the moments obtained by solving 
equations (9) and (10). 

The desired accuracy of the calculations, which is affected mainly by the choice of the discrete 
mesh, depends on the particular problem. Because the main goal of this study is the calculation 
of the spectra of accreting compact objects, the X-ray colors of which are observed to change by a 
few percent when the mass accretion rate changes by a factor of ~ 2, I require for the solutions a 
fractional accuracy of ~ 10~^ at each grid point. 

5. RESULTS 

In this section I solve the radiative transfer equation for a number of model problems with 
spherical symmetry that are related to various astrophysical systems and in particular to media 
around compact objects. In §5.1, I discuss Comptonization of soft X-ray photons by a static, uni- 
form medium of hot electrons. In §5.2. I solve the radiative transfer equation in a cold, divergence- 
free inflow; although such a configuration requires fine tuning of the physical conditions around an 
accreting object and may not occur in nature, it is used here to provide physical understanding of 
the effect of Comptonization by the bulk electron velocity. Finally, in §5.3, I discuss more realistic 
problems of spherical accretion onto and outflows from compact stars. 

5.1. Compton Scattering in a Hot, Static Medium 

For the first model problem, I consider a static, spherically symmetric, purely scattering 
medium and set the parameters to values that are typical for weakly magnetic, accreting neutron 
stars (sec Psaltis ct al. 1995). I set the outer radius of the scattering medium to 3rin, its electron 
temperature to 10 kcV, its electron scattering optical depth to 4, and the energy-dependence of 
the radiation flux at the inner boundary to be that of a blackbody with a temperature of 0.5 keV; 
the normalization of the input flux is arbitrary because I have assumed that there are no sources 
of photons in the medium. This is a frequently solved problem that has been analyzed in detail 
by, e.g., Katz (1976), Shapiro et al. (1976), and Sunyaev & Titarchuk (1980). All these authors 
showed that the emerging radiation spectrum at energies much larger than the energy of the in- 
jected photons is largely independent of the details of the input spectrum and can be described by a 
power-law with an exponential cut-off at 3Te . Here I solve this simple problem to study the fact 
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Fig. 3. — The zcroth and first moments of ttic specific intensity of tlie radiation field in the static, liot medium 
discussed in §5.1 (in arbitrary units). In botli panels the curves correspond to radii in the medium that are equidistant 
in the parameter Q (see eq.[19]). 



that the variable Eddington factors are not independent of energy even when the photon mean-free 
path is independent of energy and the medium is static (see also the discussion in Paper I). 
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Fig. 4. — The first Eddington factor at two photon energies (solid lines) as well as the ratio of the energy-integrated 
second to zeroth moments of the specific intensity of the radiation field (dashed line) for the model problem discuss 
in S5.1. 



Figure 3 shows the zeroth and first moments of the specific intensity of the radiation field, J 
and , as a function of photon energy and radius. Both quantities decrease overall with increasing 
radius because of the dilution of the radiation field. The energy dependence of the radiation energy 
density changes within one photon mcan-frcc-path from the inner boundary but remains unchanged 
at larger radii. On the other hand, the radiation flux evolves throughout the scattering medium. 
This is a consequence of the fact that the system of moment equations has been solved by imposing 
an inner boundary condition on H and an outer boundary condition on the ratio H/ J and not on 
J itself. 

Figure 4 shows the first Eddington factor /'''' = K^^ / J at two photon energies as well as the 
ratio of the energy-integrated second to zeroth moments of the specific intensity of the radiation 
field. The Eddington factor depends on photon energy even though the scattering cross section and 
hence the photon mean-free-path are mostly energy independent at these low photon energies. This 

is because photons were injected into the medium with energies <C Tg and on average gain energy 
at each scattering. As a result, photons emerging with low energy have experienced on average a 
smaller number of scatterings than photons with higher energy and therefore their distribution is 
less isotropic. 




Fig. 5. — The radiation spectrum at the outer boundary (in arbitrary units) obtained using the method of vari- 
able Eddington factors and the spectra obtained using the same boundary conditions but two energy-independent 
prescriptions of the first Eddington factor. 

Figure 5 compares the emerging radiation spectrum obtained using the method of variable 
Eddington factors with the spectra obtained using the same boundary conditions but two energy- 
independent prescriptions of the first Eddington factor that are commonly used; the second Ed- 
dington factor does not enter the calculation when the bulk velocity of the electrons is zero. Note 
that neither prescription is the correct diffusion approximation for Compton scattering because 
they both neglect the energy dependence of /'''' on photon energy shown in Figure 4. However, the 
discrepancy between the self-consistent solution and the ones obtained using prescribed Eddington 
factors is small for the prescription that depends on optical depth and can be as large as ~ 50% 
at high photon energies for the prescription that is independent of optical depth. This discrepancy 
increases when the ratio of the photon mean-free path to the smallest characteristic length scale in 
the system increases. 



5.2. Compton Scattering in a Divergence-Free Flow 



In a moving medium, the photons gain energy by scattering off the fast moving electrons. As 
measured by an observer comoving with the medium, the photon energy density increases mostly 



-16- 




Fig. 6. — The radiation spectrum at the inner boundary (dotted Hne), the emerging radiation spectrum when 
only terms that are first-order in the electron bulk velocity are taken into account (dashed line), and the emerging 
radiation spectrum that is correct to second order in the electron bulk velocity (solid line), in arbitrary units, for the 
configuration discussed in §5.2. 



due to the convergence of the flow, which is described by terms that are proportional to V ■ V 
(Blandford &; Payne 1981). On the other hand, the flux that is measured by an observer at 

infinity increases because of the systematic energy exchange between photons and electrons, which 
is described mostly by terms that are second-order in V (Psaltis & Lamb 1997). 

In order to investigate the effects of Comptonization by the bulk electron flow and distinguish 
them from the effects of the convergence of the flow, I study in this section the photon-electron 
interaction and the formation of spectra in fast, divergence-free flows. For the model problems 
shown below, I set the outer radius of the scattering medium to SOrin, the electron temperature to 
zero, and the electron velocity to the divergence-free proflle V = 0.2r~^. Finally, I set the electron 
density to a constant value, which arises from mass conservation in the flow, and characterize each 
flow by the total electron scattering optical depth in the radial direction. 

Figure 6 shows the radiation spectrum at the inner boundary, the emerging radiation spectrum 
when only terms that are first-order in the electron bulk velocity arc taken into account, as well as 
the emerging radiation spectrum that is correct to second order in the electron bulk velocity, for a 
flow with a scattering optical depth of 2. The terms that are first-order in the electron bulk velocity 
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describe mainly the effects of non-relativistic Doppler shift and hence produce a displacement in 
loge of the input spectrum towards higher photon energies (see also Psaltis & Lamb 2000). On 
the other hand, the terms that are second-order in the electron bulk velocity result in a systematic 
photon upscattering and produce a power-law tail at photon energies higher than the average 
energy of the injected photons. Note here that the approximation employed here of keeping, in 
the scattering kernel, only terms that are first order in e limits the validity of the calculation to 
energies < 100 keV. 

When the total optical depth of the scattering medium, and hence the average number of 
scatterings that each photon experiences, are increased, then on average photons gain more energy 
by scattering off moving electrons. As a result, the emerging radiation spectrum is displaced 
towards higher photon energies, and the power-law tail becomes flatter (see Fig. 7a). When the 
bulk electron velocity in the flow increases, the average energy gain per scattering increases as well, 
and hence the power-law tail also becomes flatter (Fig. 7b). As a result, scattering of photons 
by fast moving electrons produces power-law tails at high photon energies, even in divergence-free 
flows, in a way that is similar to the well understood process of thermal Comptonization. 

5.3. Compton Scattering in Inflows and Outflows 

I flnally perform simple model calculations of accretion onto and outflows from a central object 
by assuming a free-fall density profile and setting everywhere in the flow the electron temperature 
to zero and the inward radial velocity equal to the free-fall velocity from infinity onto a compact 
object with a radius of 10 km and a gravitational mass of IAMq. I also set the outer radius of 
the medium to SOrin; at these large radii the photon mean-free path is large compared to radius 
and the electron bulk velocity is small compared to the speed of light that the flow does not effect 
signiflcantly the photon distribution. Finally, instead of specifying the mass accretion rate onto the 
compact object, I specify the total electron-scattering optical depth of the flow. 

Figure 8a shows the radiation spectrum at the inner boundary, the emerging radiation spectrum 
when only terms that are first-order in the electron bulk velocity arc taken into account, as well as 
the emerging radiation spectrum that is correct to second order in the electron bulk velocity. As in 
the model calculation discussed in §5.2, the terms that are first-order in the electron bulk velocity 
describe mostly the non-relativistic Doppler shifts and produce a displacement in log e of the input 
radiation spectrum towards higher photon energies. On the other hand, the terms that are second- 
order in the electron bulk velocity systematically upscatter photons and produce a power-law tail at 
high photon energies. When the optical depth of the flow is increased, the systematic upscattering 
of photons becomes more efficient and the high-energy tail becomes flatter. 

Figure 9 shows the same information as Figure 8 but for the case of an outflow. For this model 

problems I have used the same parameters as in the inflow calculation but changed the sign of 
the velocity at all radii. In an outflow, the terms that are first-order in the electron bulk velocity 
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Fig. 7. — The dependence of the emerging radiation spectrum (in arbitrary units) on (a) the electron scattering 
optical depth and (b) the maximum electron bulk velocity, when the scattering optical depth is set to 5, for the 
configuration discussed in §5.2. 



produce a displacement in log e of the input radiation spectrum towards lower photon energies. On 
the other hand, the effect of the terms that are second-order in the electron bulk velocity is the 
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Fig. 8. — (a) The radiation spectrum at the inner boundary (dotted line), the emerging radiation spectrum when 
only terms that are first-order in the electron bulk velocity are taken into account (dashed line), and the emerging 
radiation spectrum that is correct to second order in the electron bulk velocity (solid line), in arbitrary units, for the 
inflow problem discussed in §5.3. (b) The dependence of the emerging radiation spectrum on the electron-scattering 
optical depth of the flow. 
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same as in the case of inflow, i.e., they systematicafly upscatter the photons and produce a tail 
at high photon energies. For the optical depths and bulk velocities considered here, the effect of 
the second-order terms nearly cancels the effect of the first-order terms and the major difference 
between the input and emerging spectra is the power-law tail at high photon energies. As in the 
case of the inflow, when the optical depth of the flow increases, the power-law tail becomes flatter, 
but the overall effects is signiflcantly reduced. 

6. DISCUSSION 

In the previous sections I described a method for solving numerically the radiative transfer 
equation that describes Compton scattering in spherically symmetric, static and moving media, 
and applied it to a variety of cases related to flows around compact objects. I showed that the first 
Eddington factor, which describes the degree of isotropy of the radiation field at a given point in 
space, may depend on photon energy, even though the photon mean-free-path is energy indepen- 
dent, if the photons are preferably upscattered or downscattered in energy at each scattering. I 
also demonstrated that the terms that are first-order in the electron bulk velocity alter the input 
spectrum in a way that is qualitatively and quantitatively different than the effect of the terms 
that are of second-order. 

The systematic upscattering of photons caused by bulk Comptonization is described by the 
terms that are second-order in V. This becomes apparent by the fact that, at each scattering with 
a moving electron, the average fractional energy change of a photon in the system frame is always 
positive, independent of the photon energy, and equal to 41/^/3 (see, e.g., Rybicki & Lightman 
1979, ch. 7; see also the discussion in Psaltis & Lamb 1997). If the distribution of photon residence 
times in the Comptonizing medium drops exponentially to zero at late times, then a power-law 
tail is produced at high photon energies (for total scattering optical depths > 1), in analogy to the 
generation of power-law tails by thermal Comptonization. The velocity at which the terms that are 
second-order in V become important depends on the ratio of the photon mean-free-path A^fp to 
the shortest characteristic length scale L in the problem. In general, when V{L/Xj^ip) ~ 1, i.e., in 
regions close to and within the photon-trapping radius, the terms that are first- and second-order 
in V become comparable in size (Yin Sz Miller 1995; Psaltis & Lamb 1997a) and hence the latter 
should not be neglected. 

The generalization of the method to time-independent multi-dimensional configurations is con- 
ceptually easy. However, a problem with axial symmetry requires the solution of the radiative trans- 
fer equation in five dimensions in phase space and the introduction of a large number of independent 
variable Eddington factors. Moreover, the presence of high flow velocities in multi-dimensions in- 
troduces sharp gradients in the variable Eddington factors, which cannot be guessed a priori, as 
required for the flrst iteration (Dykema, Klein, &; Castor 1996). Even with these problems, however, 
the method of variable Eddington factors will still offer a significant advantage, when compared to 
any other iterative method, which arises from the fact that the Eddington factors are bound to lie 
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Fig. 9. — Same as in Figure 8, but for the outflow problem discussed in §5.3. 



over a very narrow range of values. 

The fast convergence of the numerical algorithm presented here makes it also ideal for future 
extensions to transport problems coupled to time-dependent hydrodynamics. For the case of ac- 
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cretion onto and outflows from compact objects, the characteristic photon diffusion timescale from 
a flow of dimension R is of order 



which is signiflcantly smaller than the corresponding tdyn ~ 1 ms dynamical timescale. As a result, 

for most cases the structure of the radiation field can be obtained by solving the time-independent 
radiative transfer problem for each snapshot of the fiuid properties. Given that the variable Ed- 
dington factors will be only marginally different between successive timesteps, the convergence of 
the algorithm will be even faster than what is shown in Figure 2. 
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Appendices 



A. Moments of the Specific Intensity in Spherical Geometry 

In defining the moments of the specific intensity in the system frame I shall use the quantities 

u{z,p,e) = ^ [l+iz,p,e) +I-iz,p,e)] (Al) 

and ^ 

v{z,p,e) = -[l+{z,p,e) - r{z,p,e)] . (A2) 

The zeroth moment of the specific intensity is then 

1 r 

J{r,e) = - / u{z,p,e)dz , (A-3) 



where r depends on z and p through the relation = z'^ + p'^. Because of the assumed spherical 
symmetry, the only non-zero component of the first moment of the specific intensity is the r— 
component: 

1 r 

H'^{r,e) = ^ zv{z,p,e)dz . (A4) 



r'^ Jo 

Similarly, the non-zero components of the second, third, and fourth moments of the specific intensity 
are 





e) = 


-IT / z'^u(z,p,e)dz 
r-^ Jo 

i (J - K") , 




e) = 




e) = 
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i (iJ*^ _ Qrrr^j 




e) = 






e) = 


^ J z^u{z,p,e)dz 




e) = 


- (J - 2K'-^ + i?^'"'"'') 
8 




e) = 


I {J- 2K"- + i?'"'"'"'') 
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i?^'-^^(r,e) = i2''^'"^(r,e) = 




— {^K^^ j^rrrr^ 


B;""t><t>{r,e) = R"^'''^{r,e) = 










j^rrrr 



(A5) 



(A6) 
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Finally, I shall also use the first two Eddington factors defined by 



r = — (A7) 

^ (A8) 



B. The Generalized Source Function 



Using equation (A8)-(A12) of Paper I, I derive the generalized source function for a spherically 
symmetric system: 

S^{z,p, e) = Sf + eS^ + TeSf + VS^ + V'^Sf , 



where 



St 



S: 



± 



^4^ 



s, 



- [3 - cos2 e+(3 cos^ e-1) rn J 
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I { (-1 + ede) [3 - cos^ e+(3 cos^ ^ - l) J 
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+ (1 - ede) [5-3 cos^ ^ + (5 cos^ ^ - 3) /"^j H"^ cos 0} 
^ { [3cos2 0-1-3 (3cos2 e-1) r^] J 

+2 [1 - 3/" + (5/" - 3) cos^ d] cos 9 

+^ (-2ea, + €^5,2) [3 - cos2 e + (3cos2 B-l) r] J 

+^ (26^^ - e^dl) [5-3 cos^ + (5 cos^ O-Z) iJ"^ cos 9 
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(6 cos^ - 4) + (3 cos^ e-S)ed,-^ (3 cos^ ^ - l) e'^d'^ 
+ (2 cos^ e-l)(^l + 2ed, + ^e^a^^ 



(B6) 



where cosO = z/r and the variable Eddington factors and the moments of the specific intensity are 
defined as in Appendix A. 



C. Coefficients of the Moment Equations 

Using equations (34) and (40) of Paper I, I derive 
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C2 = -^VxSe + lvxedeSe + ^VSeedeX, (C12) 

where tq = Tr{r = rin) and I have suppressed the dependence of the various quantities on spatial 
position and photon energy. 



D. Coefficients for Differencing the Moment Equations 



Using equations (34) and (40) of Paper I I derive 
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where tq = Tr(r = r\^) and I have suppressed the dependence of the various quantities on spatial 
position and photon energy. 
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